	frames reset 
	frame create kden_est
	use "${clean_data}/panel_physicians_clean.dta", clear

	* Prepare data
	keep if year==2017
	
	foreach y of var logptotinc logaginc {
		if "`y'" == "logptotinc" 		loc label "Log-Earnings" 
		else if  "`y'" == "logaginc" 	loc label "Log AGI (1000 $)"	
	
		// Save data for histogram	
		preserve
			
			* Add percentiles 
			foreach perc in 1 5 25 50 75 95 99 {
				
				centile `y', centile(`=`perc'-0.01' `=`perc'+0.01')
			
				su `y' if inrange(`y',`r(c_1)',`r(c_2)')
				if `r(N)' < 11 {
					di "Less than 11 observations between 49.99th and 50.01st percentile for `y', use broader interval instead"	
					assert `r(N)' > 10
				}
				
				g p`perc' = `r(mean)'
				g N_UR_p`perc' = `r(N)'
			}
	
			* Add mean, sd and other locals for graph
			qui su `y', d
			g mean = `r(mean)'
			g sd = `r(sd)'
			loc p1 = `r(p1)'
			loc p99 = `r(p99)'
			loc N_UR = `r(N)'
			
			* Add data for graphv using kdensities 
			kdensity `y', nograph k(gaussian) n(200) generate(x_`y' y_`y')
			replace x_`y' = . if !inrange(x_`y', `p1', `p99')		
			replace y_`y' = . if !inrange(x_`y', `p1', `p99')		
			keep if !missing(x_`y')
			keep `y' x_`y' y_`y' mean sd p1 p5 p25 p50 p75 p95 p99 N_UR*
			g var = "`y'"
			g varlab = "`label'"
			ren x_`y' xaxis
			ren y_`y' frac
			g N_UR_kden = `N_UR'
			drbest_docinc xaxis frac mean p1 sd p5 p25 p50 p75 p95 p99
			drop `y'
			order N_UR*, last
			tempfile kden`y'
			save `kden`y''
		restore
		
		frame kden_est: append using `kden`y''
	} 
	
	frame kden_est: export delimited using "${mypath}/intermediate_csv/desc02-kden_earnings1to99_2017.csv", replace dataf delim(tab) 

   * Analysis of variance

	use "${clean_data}/panel_physicians_clean.dta", clear

	keep if year==2017
	
	eststo clear

	cap log close
	log using "${output}/descriptives/distr.txt", replace text 
	
	gen white=(race==1)

	eststo: reghdfe logptotinc, absorb(age) 
	eststo: reghdfe logptotinc female married is_foreign white, absorb(age birthplace_char)
	eststo: reghdfe logptotinc, absorb(cms_spec_code age)
	eststo: reghdfe logptotinc, absorb(cz90 age)
	eststo: reghdfe logptotinc, absorb(age einsize_discrt)
	eststo: reghdfe logptotinc with25Kbizinc, absorb(age)
	eststo: reghdfe logptotinc top5_school, absorb(age)
	eststo: reghdfe logptotinc top5_school, absorb(age cms_spec_code)
	eststo: reghdfe logptotinc with25Kbizinc, absorb(cms_spec_code cz90 age einsize_discrt)
	eststo: reghdfe logptotinc female married is_foreign white with25Kbizinc, absorb(age cms_spec_code cz90 einsize_discrt birthplace_char)
	eststo: reghdfe logptotinc female married is_foreign white with25Kbizinc top5_school, absorb(age cms_spec_code cz90 einsize_discrt birthplace_char)
	
	log close
	
	local absorb_lists ", Age, Age; Female; Married; Alien; White; POBST, Age; CMS Specialty, Age; CZ-FEs, Age; Employer size, Age; Business Inc. $>$ 25K, Age; Top5 School, Age; Top5 School; CMS Specialty, Age; CMS Specialty; CZ-FEs; Employer size; Business Inc. $>$ 25K, Age; Female; Married; Alien; White; POBST; CMS Specialty; Employer size; Business Inc. $>$ 25K, Same and Top 5 School"

	esttab * using "${output}/descriptives/desc02-earnings_variance.csv", cells(b(fmt(%09.4g)) se(fmt(%09.4g))) ///
	constant stats(N r2,fmt(0 4) labels("N_UR" "R2")) label mlabels(,depvars) ///
	eqlabels(none) ///
	postfoot("Regressors: `absorb_lists'") ///
	varlabels(_cons "Constant") replace
	
	// Add rounded N
	import delimited "${output}/descriptives/desc02-earnings_variance.csv", clear asdoub	
	insobs 1, after(18)
	
	ds v1, not
	foreach var of varl `r(varlist)' {
		loc N_UR =  subinstr(subinstr(`"`=`var'[18]'"',"=","",.),`"""',"",.)
		drbcount `N_UR' 
		replace `var' = `"="`r(rounded)'""' if `var' == ""
	}
	replace v1 = `"="N""' if v1 == ""
	
	export delimited using "${mypath}/intermediate_csv/desc02-earnings_variance.csv", replace dataf delim(tab) 
	
	
